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ABSTRACT 


Geographic location of radar emitters is the process of estimating an emitter’s 
location upon the surface of the earth from direction of arrival (DOA) data for the 
targeted emitter. The current Emitter Location (EMLOC) algorithm utilized by the 
Grumman EA-6B Prowler is based on a thesis presented ™ Mr. Richard Opperman in 
June 1982. With the advent of increased processing demands on the AN/AYK-14 
Tactical Computer as part of recent software upgrades to the AN/ALQ-99 Tactical 
Jamming System, it was hoped that a Kalman Filter, or Extended Kalman Filter based 
algorithm would reduce the processing time and memory requirements for the EMLOC 
algorithm. This thesis compares the current algorithm, and the Kalman/Extended Kalman 
Filters in a tactical scenario to determine if a change in the current Onboard Flight 


Program (OFP) should be recommended. 
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I. BACKGROUND 


A. PURPOSE 


This thesis examines current and proposed methods of geographic emitter location 
(EMLOC) with respect to the EA-6B aircraft via the AN/ALQ-99. The current algorithm 
will be compared to a proposed Kalman Filter based design. The various alternatives to 
EMLOC, along with their underlying assumptions, will be verified mathematically and a 


proposed test developed. 
B. AIRCRAFT 


The EA-6B aircraft is currently flown with three variants of the Improved 
Capabilities II (ICAP IJ) software/hardware upgrades. These are Blocks 82, 86, and 89A. 
The difference in the various blocks will be discussed, however, the focus of this paper 


will be Block 86/89A configured aircraft. 
c. AN/ALQ-99 


The AN/ALQ-99 is composed of two major sections. The Tactical Jamming 
System (TJ S) carries out the majority of the Electronic Attack (EA) capability of the 
system. The Onboard System (OBS) is effectively the Electronic Support (ES) portion of 


the system. By utilizing the superhetrodyned receivers, the operator may scan various 





prioritized portions of the electromagnetic spectrum, or the entire effective frequency 
range for the receivers. For Block 82 aircraft, as receivers are scanned, they are 
controlled by the Command Mission Computer (CMC) and the Computer Interface Unit 
(CIU). If RF energy is received by any one of the receivers, that information is sent to the 
| Encoder where information is digitized and sent to the CMC for further processing. Only 
one receiver’s data may be processed at a time. The system, therefore, has high 
sensitivity, but low probability of intercept with eight bands scanned by six receivers, and 
only one channel for processing. Block 86/89A aircraft combine the encoder and CIU to 
form a CIU/E which also assists in matching emitters with known parameters. 
Additionally, two channels are provided to the CMC to speed (in theory, double) 
processing. Test data does not support increased processing speed, however, all aircraft 
will ultimately be deployed in the Block 89A configuration and thus this version will be the 
ee of this thesis. 

As stated previously, the receivers are digitally tuned and superhetrodyned. 
Various methods are employed to prevent image frequencies from entering the system for 
further processing. However, these methods are not fool proof a operators must use 
care to employ hardware fixes, especially near high impulse density energy sources. Each 
receiver has 4-6 antennas arranged in the tail of the aircraft to localize the direction of 
arrival. For this paper we will focus iipon the six antenna arrangement. 

The six antennas are arranged at 60 degree increments to cover 360 degrees of 
space ound the aircraft. As signals are collected by the receivers, the amplitudes are 


compared at the different feeds, and a geometric average is sent to the encoder. This 





method of amplitude comparison is not very accurate, but in theory would work as 
follows: 


Assume at feed 1 the amplitude was 2 units. 
At feed 2 (clockwise from feed 1) the amplitude was | unit. 
_ Therefore, the energy is twice as great at feed 1, so the Direction of Arrival 
(DOA) is closer to feed 1, but how much closer? 
By taking a weighted average, the signal DOA must be 20 degrees to the 
right of feed 1 and 40 degrees to the left of feed 2. 
Thus if feed one is boresighted with the aircraft heading, the DOA is 20 degrees relative 


to the aircraft nose in the direction of feed 2 boresight. 


In theory, this system provides quick, accurate results, with inexpensive hardware. 
Experience, however, shows that the six antenna system achieves a nominal +/-10 degree 
accuracy. Additionally, a DOA does not guarantee sure location of the detected emitter. 
Atmospherics (refraction) and multi-path (reflections from ground obstacles) contribute 
bogus DOAs. Finally, if energy is detected at three feeds (which is possible with close 
proximity high powered signals), the detection is dropped and the receiver continues to its 


next programmed tuning step. 
D. EMLOC FUNCTION 


When a receiver detects a threat (or friendly for that matter) emitter and a match is 
made with known parameters, it is displayed to the operator in frequency and azimuth 
(after degrees true conversion). The operator may “hook” the displayed symbol and via 


several software steps, select the target emitter for Emitter Location, or EMLOC. The 


tad 





OBS then runs each subsequent detection of this emitter match through an algorithm that 
attempts to calculate the geographic location of the emitter. 

The actual mathematics necessary to accomplish this task are not terribly difficult 
to compute, but an explanation is necessary. The following chapter will describe the 
geometry and basic assumptions required to compute a target emitter’s geographic 


location given a change in aircraft position and a change in DOA. 


Il. THE PROBLEM 


A. APPLICATION 


All DOAs are delivered to the OBS relative to the aircraft. The results need only 
be developed relative to the aircraft since navigation is now precise with the advent of 
Global Positioning System (GPS) in Block 89A. Thus, with an emitter position calculated 
relative to the aircraft position (which is known with some precision) the emitter may be 
geo-located. The triangulation of emitters via Direction Finding (DF) equipment may be 


accomplished manually, and the “Flat Earth Approximation” geometry is developed quite 


simply (see Figure 1). 
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Figure 1: Flat Earth Approximation Geometry 
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tan@, = ——_—+ 2-] 
ba ees (2-1) 
and 
~X, 
tan, =—— y (2-2) 
f 
If latitude and longitude are substituted for X-Y axes, the following equations 
result: 
A A, —Ah, 
L, ~~ L, 
and 
: (A, -4,)cosL 
tan, = 4-7 (2-4) 


Pia 


where the ‘i’ subscript indicates the aircraft initial position at detection and match, and the 
‘f subscript indicates the final aircraft position. The emitters coordinates are (Xz, Y7) are 
unknown. These coordinates can be further transcribed into geographic latitude and 
longitude, but we will not pursue these calculations at this time. 

More importantly, by examining the flat earth approximation we will validate these 
“simple” equations, vice more complicated spherical solutions required by long range 
systems. Since we are only concerned with microwave frequencies, we can approximate 


the maximum distance for direct Electromagnetic propagation with: 


d= 123(Jh, + Jh, ) 


where d@ = radar horizon in nautical miles, 








h,,h, = altitude in feet for the aircraft and emitter 
assuming the 4/3 earth radius approximation for refraction characteristics. Thus if h, is 


near sea level and aircraft ceiling is less than 45,000 feet, radar horizons should not exceed 
261 nautical miles (300 statute miles). L Laddie Coburn [Ref.1] proved that at distances 
less than 360 nautical miles the planar solution was less than 0.5% in error when 
compared to the exact spherical solution, and that for ranges less than 150 nautical miles, 
the planar solution was sufficiently accurate for emitter location estimation. 

Unfortunately, the EMLOC problem is not quite this simple. With refractive 
errors and multipath, there are many arriving DOAs to consider. Also, the simple 
geometry shown in Figure 1 is not always so clear. Typical changes in bearing may only be 
2-4 degrees over several minutes flying time when aircraft aspect is nose/tail and/or range 
to target emitter is great. Therefore, various algorithms are applied that employ maximum 
likelihood estimators to locate emitters and thus, most geographic locations are a 
statistical approximation of the most probable location of the emitter. 

The algorithm currently employed by the EA-6B is an iterative method that was 
developed by Daniel Charles Opperman in June 1977. He developed his method and 
refers to it as the BOOZOO algorithm (not the term most operators would have 
preferred). During the early 1970’s, several students at the Naval Postgraduate School, 
under the tutelage of Professor Harold Titus investigated the use of Kalman Filters as a 
possible method for geo-locating emitters. We will develop both methods and evaluate 
their performance as to accuracy, processing speed and robustness under less than 


optimum geometry. 





I. THE BOOZOO ALGORITHM 


A. THE BOOZOO ALGORITHM DEVELOPED 


As stated previously, the Boozoo Algorithm (BZA) was developed by Opperman 
as part of his thesis work with the Tactical Electronic Reporting, Processing, and 
Evaluation System (TERPES) [Ref. 2]. This system was developed by the U.S. Marine 
Corps to take ES data collected by a variety of sources (primarily the EA-6B) and develop 
the Electronic Order of Battle (EOB). The system is currently fielded and is organic to 
the Marine Tactical Electronic Warfare Squadrons. Opperman developed his thesis based 
upon equations provided by the North American Rockwell Corporation (NAR). 
Opperman noted that the NAR equations were based on a model where the ES aircraft 
flew a track long enough to create at least 90 degrees of change in bearing. This is not 
always possible so the BZA was developed to account for DOA oe that was significant 
when compared to the change in the initial bearing and final bearing. The BZA was 
eventually incorporated into the TERPES system in the mid 1980’s and into the EA-6B 
aircraft during the early 1990’s with the EMLOC function. 

Opperman used the following as a foundation: 


P(q@, \da, = oe), ) da, (3-1) 


where P(a@,) is the probability that the ith aircraft position will register a bearing between 


6,+a, and 6, +a, +da,, and a,is the DOA measurement error for the ith 








measurement. The value, o, is the standard deviation of the DOA bearing error 
associated with the /th andi 

If we take n independent measurements, then we have joint probability: 
exp(-1/2 f(E)) 


P(a,,...,@,)da,da,...da, = a 
(27) ‘\o,0,...6, 


da,d,...da, (3-2) 


where 


s(E)=> ws ) on 


/ 

and (see figure 2) E is the most likely position of the target emitter, 7 is the aircraft 
position vector, and ¢ is a vector miles to 7, and lying in the plane defined by the 
vectors, Q, (actual position vector) and E , 


E Plane 









Earth Surface (Great Circle) 


Possible Choice for 
Position of Target Emitter 


Cutting Plane 


Figure 2: Unit Sphere Approximation of Emitter Location 
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In order to minimize the error, Opperman suggested a weighting function that 
would be changed according to where the DOA came from. In our case this would imply 
some DOAs are more accurate than others. This turns out to be true in actuality and 


therefore a weighting function is defined: 


| 
We = Veg? (3-4) 


where & 1s an arbitrary constant of proportionality. 


Thus the probability would be maximized if we minimize 
f(E) =k Wa? (3-5) 
i=] 


As part of his foundation, Opperman does this, by first making an approximation. If @, is 


small then sing, ~ @, So, 


f(E)= ky W? sin? a, (3-6) 


i=] 


From Figure 2, 


sing, ~ col - 0, | / sin B, 


so now, 
f(E) = 2 W? (cos* (7% - o,)/sin® B,] (3-7) 
We may assume that 
. 9 lo. > 
sin® B, x —> sin’ B, 
A jz 


and making another substitution with simplification, 








n 


; > Ww? cos"(7/ a] 

f(E) = nk #1 _____— (3-8) 
| > sin’ B, 

i=] 


Equation (3-8) when viewed with Figure 2 may be expressed in terms of the 


aA 


expected emitter position, EF . 


f(E)=n anes (3-9) 








The function in equation (3-9) may be minimized, but some other reference may 
need to be defined. Assuming / and ¢ are longitude and latitude of the aircraft position 
respectively, and @ remains the measured DOA as measured from true north, we may 
redefine the vectors, Ef, and ¢ such that: 


a, = cosa, cos¢, 
6b, = sind, cos¢, 


c, =sing, 
X, = sind, cos@, —sing, cosA, sing, 
Y, = —cos4,cos6, —sing, sind, sind, 


Z, = cos¢@, sind, 


Where. Finally, we use earth centered coordinates (i,j,k) to define: 


Now substituting into equation (3-9) 


f(E)= gh (3-10) 
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where 


g= nk W(X +¥v+Z,w) 


i=] 
h= > |e —¢v) +(cu-a,w) +(a,v- b,x) | 
i=] 
By differentiating with respect to u, v, and w we obtain the gradient of the function. 


Vf ( tf (Z :)) will have three components 


f, =(g, ne v, Gi} 


Now, if the gradient is allowed to approach zero, then f (E) will be minimized 


and equations (3-11) will become 


g, =h,f(E) 
g, =h,f(E) (3-12) 
g. =h,f(E) 


By differentiating functions g and / and substituting into equations (3-12) we 


obtain: 
Qi, @. Az] U b, 5, 3 | u 
@, 4x2 4g3 || V = T(E) by, By. ba3 | V (3-13) 
a3, 43, 33 |W b,, 63, 63; | w 

where 





o3.> DAW SE aes DAZ W, 
i=} i=) 
a3 = > X,Z,7 a3, = Af 


and 


The solution for u, v, and w is simplified if we assume that the right side of 


equation (3-13) has a known value. Therefore by using a method of iteration we obtain: 


GQ, @. 3 | Uy, b, 5, 5, | u, 
Qy Ay Ay | Via) | = r(E) Dis Dy. De Vy, (3-14) 
Qi3 Ap3 33 | Wy) b; 5; 63; | w, 


where k = 0,1,2,3,...m, and the zeroth iteration of wu, v, and w are the solutions for the 


emitters position. 


Simplifying the right half of equation (3-14), 


Uys) Q, Ay 3 |W 4, B, 
Al Vy) |=]. Ay 5 | Yea = f(E} B, (3-15) 
Wray G3 4.3 433 | Wy, B, 


where 


B, = 6,,u, +b,.v, +b,,w, 
B, =b,,u, +bv, +5,,w, 
B, — b, ,u, + b,,v, + b,,W, 
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Equation (3-15) can be solved by Cramer’s Rule, provided the A matrix is singular. 











Then, 
f(E} B, Qi. Ay; 
Uy = ree det|B, a, 4), 
: B; aj; ay, 
f(E) a, B, ay; 
Vigo ry detia,, B, a,, | (3-16) 





4 : detla,, a, B, 


This solution, however, leads to the determination of the position vector in both 


direction and magnitude. Since we only need direction, we may further simplify the 


solution. 
By 4, 4; 
U,,, =| det}B, a, a, 
By, a; ay; 
a, B, a, 
Ven =| deta, B, a, (3-17) 


The remaining calculations are quite simple to find the latitude and longitude of the 


A 


emitter position vector, E . 
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Longitude = arctan(v,, /1,, ) 


Latitude = arta" | 
2 2 
V ui, a se 


The BZA appears to be a math intensive algorithm that would require an 
inordinate amount of processing time when compared to other methods of estimation. It 
is NOW appropriate to develop several alternative estimators to the same problem, the 


Kalman Filter, and Extended Kalman Filter. 
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IV. THE ALTERNATIVES 


A. THE KALMAN FILTER 


The Kalman Filter equations are a mathematical representation of linear systems 
and are adaptive to a variety of situations when trying to simulate, or estimate linear 
systems. They have extension to nonlinear systems, but may not be as accurate in certain 
circumstances. The application of Kalman Filter equations to the EMLOC problem have 
been explored since the early 1970’s. In June, 1972, L. Laddie Coburn [Ref. 1] developed 
two algorithms, a Kalman Filter and an Extended Kalman Filter, since there are 
nonlinearities. Also of use was the follow-on work to Coburn by Edward H. Mills [Ref. 
3] that was useful in creating algorithms compatible with the AYK-14. We will examine 
both algorithms. 

The underlying mathematical model for the Kalman Filter tracker consists of a 
State equation (a linear difference equation) and a linear measurement equation. Noise 
processes appear in both equations. 

The state equation is: | 

X(k + 1) = O(k + Wk) X(k) +TW(k) (4-1) 
where ®(k + I|k) is the state transition matrix, T is a distribution matrix related to the 
random forcing function, W(k) is a random forcing function that accounts for random 
excitations (noise), and k is the time step. 


The measurement equation is: 


Z(k) = H(k)X(k)+V(k) k= 1,2,3,... (4-2) 
where H(k) is the observation matrix, and V(x) is the measurement noise. 
The noise is assumed to have zero mean and therefore, 
E[V (kV Ca)" | = REDO) 
TE[W(R)W IC = O(b)6(¥) (4-3) 
E[V Wa)" | = 0, Vk, j 


where 0(kj) is the Kronecker delta function and has value zero unless k =. 


The Kalman recursion equations are summarized: 


G(k) = PCklk —1)H(k)"[H(k) P(k|k — DAK)! + RGD] 

P(k\k) = P(k|k - 1) - G(K)A(k) P(Alk -1) 

P(k + 1k) = O(K +1, k)P(AK)O(K +14)? +O(k) (4-4) 
X(kik) = X(k\k —1) 4 G(k)[Z(k) ~ H(k) X (kk - | 

X(k{k -1) = O(k,k -1)X(k - Wk -1) 


where X (| 7) is the estimate of state X(k) based on j measurements (Z(1), Z(2),...Z(j)). 
The Kalman filter gains are represented by matrix G(x) and P(k|/) represents the error 
covariance matrix of the estimates. 

Since we are concerned with DOAs and bearing rates, we will filter noisy observed 


DOAs to obtain an optimal estimate. This is best approximated by a state transition 


matrix: 


| 1 7(k+}) 
O14) =| | 


(4-5) 
H(k)=[1 0] 





Also, since DOAs are the only observed value, R(k) and W(k) become scalars and the 
bracketed terms in equations (4-4) become scalars, also. This will play a part in 
simplifying and reducing the matrix equations to scalar form. The non-uniformity of 
arrival time for DOAs requires that the term 7(k+1) remain a variable in each of the 
applicable recursion equations. This requires that the error covariance and gains be 
computed continuously since they are a function of the sample interval. 


The P and Q matrices have scalar components given by 


ey 2 O 
re a “al and the G matrix is a2 x 1 matrix, bs 
Py Py QO, 0». G, 


Rewriting the gain equation G(A), 


G)_[P FeP ly, op’ TT, pl - 
ol FP, Py 0 | iB elo ; 


The observation matrix H(k) = [1 0] forces the inverse term to become a scalar, allowing 


us to compute gain terms directly: 


Gi} Fy -1 
ehh ° 


which results in scalar gain equations 


a P,, (Alk -1) 
COP Elk 1) + RO) a 
G, (k) = P,,(k\k — 1) 


P(e -1) + R(k) 
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Now the first DOA is not filtered, but a technique of filtering backward is available 
to improve the estimate of the initial DOA. These equations are similar to the Kalman 
equations, but they update and smooth previous data based on the last DOA that has been 


optimally estimated. This method was proposed by Rauch [Ref. 4] in 1963. 


X (ik) = X(k -1) + D(WA)G(K)] Z(H) ~ A(k)X (kk - 1)| 


(4-9) 
D(|k) = DQ|k - 1) P(A - Wk — 1) O(k, k - 1)" PRK 1)" 
where D(1|k) becomes a function of 7(k) and must be solved sequentially before the 
6(k) = H(k)O(k) +V(k) (4-10) 


where D(1|k) becomes a function of T(x) and must be solved sequentially before the 
smoothed estimate can be computed. These equations can also be reduced to scalar form 


and were utilized by Coburn in his thesis work. 


Initialization of the filter requires knowledge of the system dynamics, statistical 


properties of the filter, and information regarding the initial state of the system. 


The noisy observation of the state may be made by 


Sie bea 


Ok) 
rtd 


= bes 


where 6(k) is a vector of noisy observations, and O(k) is a state vector of exact emitter 


bearing angle and bearing rate. The observation matrix H(k) = [1 0] is for measurement of 


the bearing angle only. Bearing rate depends on speed, range and heading of the aircraft 
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with respect to the emitter. Since the range is initially unknown, it must be estimated. 
This could be a standardized variable based on the emitter to be located, or could be an 


operator input. For this exercise it will be denoted by the variable r. Therefore, 


7.29578) / | 
a rad 
= - sin(Relative DOA) a ae deg/ sec (4-11) 


The remaining equations are simple substitutions of theta for X in the basic Kalman 


equations presented, resulting in the optimal estimation equation: 


(mk) ]_[1 THY Ok-1k-Y] (GH) 
ace) : E lie 1k - 7 . CO eE® (4-12) 


where 


E(k) = (k) - 6(k ~ Ik -1) — T(A)O(k - 1k - 1) 


Again, since there is not a uniform sampling interval, the filter gains G may vary 
and not approach a steady state value uniformly. 
Once the values for theta are determined, an iterative approach to the equations in 


chapter two should lead to the emitter position. 
B. THE EXTENDED KALMAN FILTER 


Due to the nonlinearities of the trigonometric functions, theoretically an Extended 


Kalman Filter (EKF) is warranted for our situation. Realistically, we will not have a 


21 





priorii em rate, nor will our initial guesses at bearing rate be accurate. Therefore, we 
must face reality and attempt to locate with only DOA data. If this is the case, however, a 
nonlinear transformation must be utilized to transform observable DOAs into filtered 
position estimates. 

This is not as intimidating as it may sound, using the Extended Kalman Filter 
(EKF) with the help of small perturbation theory and a Taylor Series about an initial point. 
Equation (4-1) then becomes 

Z(k) = N[X(k)] + V(A) (4-13) 

where Z(k) represents observable DOAs, X(k) is the emitter position vector, and NV 


represents the nonlinear transformation 


és A,—A 1B 
@ = rctan{ x —AJeoste (4- 14) 
L,-L 
If the position error of the state vector X(k) is given by 
X(k) = X(k)— X(k), (4-15) 


then the true position X(k) may be expanded about the most recent optimal estimate, 


X (A) in a Taylor series expansion with higher order terms thrown away as shown: 


N[X(k)] = N[X(A)] + M(A) XH... (4-16) 
where 
aN j ; 
Mw)= 3). =|Y%a, Ya, 





for the emitter location algorithm given in equations (2-3) and (2-4). 


The Kalman filter recursion equations become 
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G(k) = PCR|k -1)M(k)"[M(k)P(R|k - Me)? + ROD] 
P(k|k) = P(k|k -1) - G(k)M(k) P(k\k -1) (4-17) 
P(kK +1)|k) = O(K +1,k)P(R|K)O(k +1,k)7 + O(K) 


Since the emitter is almost always stationary (or near stationary), the ® matrix becomes 
the identity matrix, greatly simplifying the equations. 

EKF filters can be volatile, and without proper initialization will diverge. 
Obviously, if the equations can not be initialized a majority of the time with good results, 
the EKF will be of no use to the operator. 

Now that the various algorithms, both current and proposed, have been described 
mathematically, it is time é subject each to the same scenarios and determine which, if 


any, is clearly more reliable, or efficient. 
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V. THE SIMULATION 


A. TESTED PARAMETERS 


The ultimate goal of this project is to program the three algorithms presented to 
determine which method: 

1.) is more computationally efficient (fewest floating point operations, or flops) 

2.) requires least amount of memory (smallest matrices and arrays held in RAM) 

3.) converges to an acceptable range of the target emitter in the fewest number of 


iterations. 
B. SCENARIO 


A realistic simulation was developed to test the algorithms using Monte Carlo 

| techniques. All computer code was written in MATLAB. 

The simulated ES aircraft flies North along the Prime Meridian (0 deg Longitude) 
from 1.5 deg S. Latitude to 1.5 deg N. Latitude, a distance of 180 Nautical Miles (Nmi). 
The aircraft will flies at constant altitude and heading at a true airspeed of 360 knots. 
Measurements are taken at 12 second intervals in an attempt to geo-locate the target 
emitter. The true emitter location is on the equator (0 deg latitude) at 1.5 deg East 
Longitude Direction of arrival accuracy was initially assumed to be +/- 10 degrees. A 


representation of each algorithm is included in the Appendices. 
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Also, in an effort to determine if geometry and/or sampling interval has an effect 
on any algorithm’s performance, the start and stop points were altered. 
The algorithms were tested as described above, and a chapter is devoted to the 


results of each. 
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VI. BOOZOO RESULTS 


A. BOOZOO ALGORITHM 


The conversion of Opperman’s code was completed with the assistance of Mr. 
Lester Hathaway, NAWC Weapons Division, NAS Pt. Mugu, CA. Additionally, the 
actual weighting scheme for one of the frequency bands was employed. Figure 3 is a 


graphic representation of the weights applied to each arriving signal. Although every 


‘ 
AY 


L 0 
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Figure 3: EA-6B Mid-Frequency Weighting Scheme 
effort was made to economize the volume of computer code written in support of the 


algorithm, the AN/AYK-14 does not support matrix operations and must therefore rely 
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upon scalar equations. As performed for this simulation, the BZA requires memory 
positions for 75 variables, including 37 weights from the weight look-up table. An 


example of the code utilized is included as Appendix A. 


B. BASIC SCENARIO 


Figure 4, depicts a typical ES track along the Prime Meridian and the cluster of 


estimated positions for the target emitter for one run of the simulation. Accuracy was 


Latitude (deg) 
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Figure 4: Single Run of BOOZOO Baseline Simulation 
impressive as shown in Figure 5 for both one simulation run, and the average of 1000 
simulations. The number of iterative loops per DOA cut was constant or nearly constant 
at two iterative loops with each loop comprising 42 floating point operations. Set up for 


the iterative process took an additional 68 floating point operations for a total of 152 
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floating point operations per typical DOA cut. Finally, a simple time analysis was 
accomplished by evaluating the time to complete one thousand iterations of the complete 
scenario. The one thousand simulations took 2,256.2 seconds to complete, for an average 
of 2.5262 seconds for each loop through the sosiitio: Thus, with 150 measurements 
taken, the BZA process averages 0.0166 seconds per measurement to arrive at a position 
estimate. Figure 5 depicts the average error associated with the BZA under the conditions 
set in Chapter V. Figure 6 (shown on the next page) is a blown up depiction of average 


error focusing upon the area of convergence for the targeted emitter. 
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Figure 5: Averaged Error 1000 Simulations 


The downward slope of the average error is very steep starting at about 5O measurements, 


and flattens out after 75 measurements, which corresponds with the aircraft passing the 
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abeam position to the target emitter. Measurements were terminated at 150, which 
corresponds to a theoretical 90 degree change in target position. Final average range error 


is less than four nautical miles (Nmi.) at the termination of the test. 
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Figure 6: Averaged Error (<20 Nmi) 1000 Simulations 


C. EFFECTS OF SAMPLING INTERVAL 


A more realistic measurement interval in a “tactical” situation would be on the 
order of 30-60 seconds between measurements. Therefore, the BZA was tested as in the 
previous scenario and average errors are depicted in Figure 7 and Figure 8 (page 31) for 
30 second interval and Figure 9 and Figure 10 (page 32) for the 60 second interval. In 


both cases, the steep slope is evident at about one third of total measurements with 
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Figure 7: Averaged Error 30 Second Intervals 1000 Loops 
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Figure 8: Averaged Error(<20 Nmi) 30 Second Intervals 1000 Loops 
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Figure 9: Averaged Error 60 Second Intervals 1000 Loops 
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Figure 10: Averaged Error (<20 Nmi) 60 Second Intervals 1000 Loops 
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a leveling off following the passage of the abeam position. Accuracy at the completion of 
the scenario was under 6 Nmi. for the 30 second interval and about 7 Nmi. for the 60 
second interval. Obviously, more data results in more accuracy, but in the three cases 
(Basic, 30 sec., 60 sec.) presented, best accuracy is achieved just after 45 degrees of 


relative change in DOA. 


D. EFFECTS OF START/STOP POSITION 


In an effort to determine if the algorithm was affected by changes in position for 
start and stop of the measurements, the basic scenario was modified so that sampling 
interval remained 12 seconds, but the initial and end points were modified. Figures 11 
through 13 (pages 34-35) show the effects of starting at 2 degrees S. Lat. And continuing 
North to the Equator for a total of 120 measurements. .Figures 14 through 16 (pages 35- 
36) depict the effects of starting at the Equator and continuing North to 2 degrees N. Lat., 
also for 120 measurements. The algorithm gave similar final results in both cases, for an 
average error at the conclusion of the scenario of about 6-7 Nmi. There is, however, a 
striking difference in the slope of the range error curves, the range error for the Equator 
start is much steeper and approaches near steady state by the thirtieth measurement. In 
contrast, the Equator finish scenario (Figure 12) requires over 60 measurements to 
achieve the same average accuracy. Clearly the preferred method to locate the target if 
the 90 degree option is not available is to start with the target at least abeam the aircraft 
and to continue on a course perpendicular to the radius of an arc described by the 


theoretical range to the target. This is due in part to the more rapid change in apparent 
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Figure 11: 2 Deg S. - 0 Deg Single Run Estimates 
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Figure 12: 2 Deg S - 0 Deg Averaged Error 1000 Loops 
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Figure 13: 2 Deg S - 0 Deg Averaged Error (<40 NMi)1000 Loops 
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Figure 14: 0 Deg - 2 Deg N. Single Run Estimates 
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Figure 15: 0 Deg - 2 Deg N. Averaged Error 1000 Loops 
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Figure 16: 0 Deg - 2 Deg N. Averaged Error (<20 Nmi) 1000 Loops 
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bearing from the aircraft to the target when the target is at the abeam position. Figures 11 
and 14 show the “memory” of the algorithm where it becomes apparent that the algorithm 
relies heavily upon the first DOA as all estimates seem to come at or near the first relative 


bearing line. 
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VIT. KALMAN FILTER RESULTS 


A. KALMAN FILTER 


The Kalman Filters tested for this thesis are representations of those developed 
and already cited by Mills and Coburn. Coburn’s thesis relied on an inverse matrix 
manipulation that is not supported by the AYK-14. His example is presented as a 
comparison to the actual tested example, as presented by Mills. Mills developed a scalar 
approximation to the inverse of the P matrix. The same baseline scenario was run with 


both algorithms an yielded the results presented in Figure 17. 


Mills Katman Filter 
oO © Coburn Kalman Filter 





Figure 17: Comparison Between Mills and Coburn Kalman Filter 
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Thus, although Coburn’s case is a clear winner, we cannot use his example since it cannot 


be supported by the current computer utilized by the typical ES aircraft. 


B. BASIC SCENARIO 


Figure 18, depicts a typical ES track along the Prime Meridian and the cluster of 


estimated positions for the target emitter for one run of the simulation. Accuracy was 
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Figure 18: Single Run Kalman Filter Baseline Simulation. 


not as impressive as the BZA and is shown in Figure 19 for both one simulation run, and 


the average of 1000 simulations. The algorithm utilizes 115 floating point operations per 
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measurement which is almost a 25 percent reduction in necessary calculations when 
compared to BZA. Finally, a simple time analysis was accomplished by evaluating the 
time to complete one thousand iterations of the complete scenario. The one thousand 
simulations took 1657.2 seconds to complete for an average of 1.6572 seconds for each 
loop through the scenario. Thus, with 150 measurements taken, the Kalman Filter (KF) 
process averages 0.011 seconds per measurement to arrive at a position estimate, which is 
almost 34 percent faster than the BZA. Figure 19 depicts the average error associated 
with the BZA under the conditions set in Chapter V. Figure 20 (shown on the next page) 
is a blown up depiction of average error focusing upon the area of convergence for the 


targeted emitter. 
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Figure 19: Averaged Error 1000 Simulations 
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The downward slope of the average error is very steep, starting at about 50 measurements 
and flattens out after 75 measurements, which corresponds with the aircraft passing the 
abeam position to the target emitter. Measurements were terminated at 150, which 


corresponds to a theoretical 90 degree change in target position. Final average range error 


is around 25 Nmi. at the termination of the test. 
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Figure 20: Averaged Error (<60 Nmi) 1000 Simulations 


C. EFFECTS OF SAMPLING INTERVAL 


As in Chapter VI, a more realistic measurement interval in a “tactical” situation 


would be on the order of 30-60 seconds between measurements. Therefore, the Kalman 
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Filter was tested as in the previous scenario and average errors are depicted in Figure 21 
and Figure 22 (page 44) for 30 second interval and Figure 23 (page 44) and Figure 24 
(page 45) for the 60 second interval. In both cases, the steep slope is evident at about one 
third of total measurements with a leveling off following the passage of the abeam 
position. There is a spike in both scenarios at around the abeam position, but in both 
cases the algorithm rapidly returns toward convergence on a minimum error. Accuracy at 
the completion of the scenario was about 26 Nmi. for the 30 second interval and about 28 
Nmi. for the 60 second interval. In the case of the Kalman Filter it appears, more data 
results in more accuracy, but not a remarkable accaeniaiai when compared to the 
baseline. Again, in all three cases (Basic, 30 sec., 60 sec.) presented, best accuracy is 


achieved just after 45 degrees of relative change in DOA. 
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Figure 21: Averaged Error 30 Second Intervals 1000 Loops 
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Figure 22: Averaged Error (<60 Nmi) 30 Second Intervals 1000 Loops 
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Figure 23: Averaged Error 60 Second Intervals 1000 Loops 
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Figure 24: Averaged Error (<80 Nmi) 60 Second Intervals 1000 Loops 


D. EFFECTS OF START/STOP POSITION 


In an effort to determine if the Kalman Filter algorithm is affected by changes in 
position for start and stop of the measurements, the basic scenario was modified so that 
sampling interval remained 12 seconds, but the initial and end points were modified. 
Figures 25 through 27 (pages 45-46) show the effects of starting at 2 degrees S. Lat. And 
continuing North to the Equator for a total of 100 measurements. Figures 28 through 30 
(pages 47-48) depict the effects of starting at the Equator and continuing North to 2 
degrees N. Lat., also for 100 measurements. Both algorithms gave similar final results for 


an average error at the conclusion of the scenario of about 35 Nmi. There is, however, a 
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striking difference in the slope of the range error curves, the range error for the Equator 
start is much steeper and approaches near steady state by the fiftieth measurement. In 
contrast, the Equator finish scenario (Figure 26) requires over 80 measurements to 
achieve the same average accuracy. Clearly the preferred method to locate the target, if 
the 90 degree option is not available, is to start with the target abeam the aircraft and to 
continue on a course perpendicular to the radius of an arc described by the theoretical 


range to the target. This is in agreement with results obtained utilizing BZA. 
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Figure 25: 2 Deg S. - 0 Deg Single Run Estimates 
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Figure 26: 2 Deg S - 0 Deg Averaged Error 1000 Loops 
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Figure 27: 2 Deg S - 0 Deg Averaged Error (<100 NMi)1000 Loops 
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Figure 28: 0 Deg - 2 Deg N. Single Run Estimates 
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Figure 29: 0 Deg - 2 Deg N. Averaged Error 1000 Loops 
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Figure 30: 0 Deg - 2 Deg N. Averaged Error (<60 Nmi) 1000 Loops 


Figures 25 and 28 show the “memory” of the algorithm where it becomes apparent that 
the algorithm relies heavily upon the first DOA as all estimates seem to come at or near 


the first relative bearing line. This is expected with the Kalman Filter as the first DOA is 


continuously smoothed as the algorithm proceeds, and is consistent with BZA. 
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VIN. EXTENDED KALMAN FILTER RESULTS 


A. THE EXTENDED KALMAN FILTER 


The extended Kalman Filter was programmed to mirror the work of Coburn and 
Mills, who both relied upon a more accurate measuring device then the one simulated for 
this scenario. Analysis of the algorithm was incomplete due to the volatility of the 


results. 


B. BASELINE SCENARIO 


The results of the baseline scenario were so poor that an attempt was made to see 
if the problem was numerically unstable due to the use of latitude and longitude. The use 
of the X-Y plane could facilitate solving the equations and a conversion based on aircraft 
present position and heading could be used to arrive at a solution. Unfortunately, for DF 
accuracy greater than 1 degree, the filter always diverged. An attempt was made to 
provide only exact initialization data to start the algorithm, with like results. Thus, there 


is no apparent utility in this method at this time. 
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IX. CONCLUSIONS 


Although the stated objective of this thesis was to determine if the Kalman Filter 
algorithms would be faster and more efficient than the current implementation, it is noted 
that accuracy is the ultimate goal and should not be sacrificed in order to increase speed, 


and efficiency. 
A. EFFICIENCY 


Clearly, the evidence presented in chapters VII and VIII shows that the Kalman 
Filters are faster and require fewer floating point operations. The AN/AYK-14 is a time 
share computer that has to be optimally programmed to expedite all tasks; however, | 


accuracy remains paramount. 
B. MEMORY REQUIREMENTS 


Again, the evidence shows the Kalman Filter requires almost 25 percent fewer 
memory registers in order to complete the necessary equations to effectively locate target 
emitters. On the other hand, the number of required memory registers involved are 
insubstantial, and the case can be made that this is opal important criteria for 


determining which algorithm is best utilized. 
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C. ACCURACY 


Without doubt, the deciding factor for aircraft operators given the limited 
disparity in time and memory requirements is effective accuracy. The BZA is clearly far 
more accurate than any of the Kalman Filters. Thus, despite the increased processing 
demands placed upon the AN/AYK-14, the BZA should be retained. The BZA was more 
accurate in virtually every scenario encountered. The effects of geometry and sample 
interval were noticeble, but were still better than both Kalman Filters (Coburn and Mills) 


in the baseline scenario. 


D. FUTURE DEVELOPMENTS 


As the algorithms present themselves now, the BZA is the clear winner. With the 
advent of portable terminals that interface with aircraft systems, the lack of speed and 
computing power now experienced will be overcome by the processing power of 
microchips. It appears that for this situation at this time, the Kalman Filter estimators are 
not the answer. That is not to say a more complicated Kalman Filter used in conjunction 
with an Extended Kalman Filter could not be developed, but the BZA is more accurate 


and not remarkably slower than the currently programmed examples evaluated. 
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E. AREAS FOR FUTURE RESEARCH 


As part of the comparison of the various algorithms, the scenarios and variables 
have been kept very simple, but true inflight acta has yielded results that are almost 
never as accurate as those presented in this paper. I would propose that the effect of 
utilizing the BZA with an 8-bit computing system inflight vice the 32-bit system utilized 
for this thesis be investigated. Also, I have kept the sample interval regular and not 
extreme in length. If real world constraints and limitations are considered, the results 
may see some impact. Additionally, I have not considered the various time averaging 
program steps the ALQ-99 applies to incoming DOAs when the time between reception is 
long. 

The EA-6B is the planned C2W aircraft for the next twenty years, and as such, 
every effort must be made to optimize the ES capability of the crew as the threat advances 


in technology and processing power. 


55 











. 36 











APPENDIX A. BOOZOO CODE 


This appendix is a representative case of Opperman’s code FORTRAN written in 
MATLAB. 


%oTo Vo Vo % BOOZOO.M %o%o To To% 

To —  %o% 
% % 
% Maj S.P. Jones, USMC 

% ‘This script is a MATLAB adaptation of the FORTRAN subroutine 

% program developed by Daniel Opperman in June, 1997. 


% flong is aircraft position longitude 

% flat is aircraft position latitude 

% azz 1s the raw doa azimuth bearing measurement 

%o wt is the weighting factor for the measurement 

% 11s the index denoting the number of the present measurement 
% ncuts is the total number of measurements 

% elong is the estimated emitter longitude in radians 

% __ elat is the estimated emitter latitude in radians 

clf 

%%% Initiate the program 

randn(‘seed',2) % set seed generator 

pirad = pi/180; % conversion factor from deg-->rad 
flong = 0; % 0 degrees long (rad) 

flat = -1.5*pirad; % 1.5 deg S. (rad) 

azzout = []; 

flatout = []; 

flongout = []; 

elatmat = []; 

elongmat = []; 

itout = []; 


tgtlatmat = []; 
tgtlongmat = []; 
range_error = []; 


tgtlong = 1.5*pirad; % 1.5 deg E. (rad) 
tgtlat = 0; | % at the equator 
wt=1; 


wttable = [.8 .96 .87 .74 .51 .5 .47 .46 .48 .57 .57 ... 
.62 .42 .64 .61 .64 .62 .67 .67 .65 .6 .62... 
.63 .7 .39 .58 .51 .45 .47 .47 57 .78 .88 ... 
.88 .99 .78 .47]; 

index = [0:149]; 


for i= 1:150 


dferror = 10*pirad*(randn); 
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Yoazz\ = atan2((tgtlat-flat),(tgtlong-flong)) 
azz = (pi/2)-atan((tgtlat-flat)/(tgtlong-flong)) + dferror; 
azzout = [azzout,azz]; 


%%% DETERMINE WEIGHT 
if azz < 180 
wtindex = round((azz + 2.5)/5); 
else 
wtindex = round((360 - azz + 2.5)/5); 
end —%%% if azz<180 
wt = wttable(wtindex); 


if i==] 
test = 2; 
% 
%o% INITIALIZE A AND B MATRICES IF THIS IS FIRST MEASUREMENT%% 


all =0; 
al2 =0; 
al3 = 0; 
a22 = 0; 
a23 = 0; 
a33 = 0; 


bli =0; 
b12 = 0; 
b13 = 0; 
b22 = 0; 
b23 = 0; 
b33 = 0; 


end %end if i== 
%% CALCULATE DIRECTIONAL NUMBERS XS,YS,ZS for Ith CUTTING PLANE %% 


slat = sin(flat); 

clat = cos(flat); 
slong = sin(flong); 
clong = cos(flong); 
$aZz = sin(azz); 
Cazz = cos(azz); 


XS = (slong*cazz - slat*clong*sazz) * wt; 
YS = (-clong*cazz - slat*slong*sazz) * wt; 
ZS = clat*sazz*wt; 


J%oTo%o CALCULATE THE A MATRIX 


all =all + XS42; 
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al2 =al2 + XS*YS; 
al3 = al3 + XS*ZS; 
a22 = a22 + YS‘2; 
a23 = a23 + YS*ZS; 
a33 = a33 + ZS42; 


% Determine the vector components X,Y,Z of A/C position vector 


X = clong*clat; 
Y = slong *clat; 
Z = slat; 


% Calculate the B matrix 


bli =bl1 + Y42 + Z/2; 
b12 = b12 - X*Y; 
b13 = b13 - X*Z; 
b22 = b22 + XA2 + ZA2; 
b23 = b23 - Y*Z; 
b33 = b33 + X42 + Y42; 


%%%% Calculate the Co-Factors of the A Matrix 


cll = a22*a33 - a23“2; 
c12 = a13*a23 - al2*a33; 
c13 = a12*a23 - al3*a22; 
c22 = al1*a33 - al3%2; 
c23 = al2*al3 - al1*a23; 
c33 = al1*a22 - al242; 


u=cll; 
V¥=clzZ: 


w=cl3; 


elong = 999; 
elat = 999; 


ifi<2 %%% Need at least two cuts 
elat = 0; 
elong = 0; 

end % If i<2 

iterations = 0; 

if i == 150 
flops(0) 

end 


while test > 1.00000000 %%%lterative Loop 


iterations = iterations +1; 


elon = elong; 
elat1 = elat; 
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Al =b11*u + b12*v + b13*w; 
A2 = b12*u + 622*v + b23*w; 
A3 = b13*u + b23*v + b33*w; 
%%% DETERMINE COMPONENTS OF POSITION VECTOR 
u=Al*cll + A2*cl2 + A3*cl13; 
V=A1*cl2 + A2*c22 + A3*c23; 
w = Al*cl3 + A2*c23 + A3*c33; 
%%o% DETERMINE LAT AND LONG OF TGT EMITTER 
elong = atan(v/u); 


elat = atan(w/sqrt(u*2 + v2)); | 
test = abs(elonl-elong) + abs(elat1-elat); 


end %%% (while test >1e-8) 

if 1 ==150 | 

flops 

end 

itout = [itout, iterations]; 

test = 2; 

flatout = [flatout,flat]; 

flongout = [flongout,0}; 

flat = flat + ((6*12/60* 1/60)*pirad); 


tgtlatmat = [tgtlatmat,tgtlat/pirad]; 
tgtlongmat = (tgtlongmat,tgtlong/pirad]; 


elatout = elat/pirad; 
elatmat = [elatmat,elatout]; 


elongout = elong/pirad; 
elongmat = [elongmat,elongout]; 


error = sqrt((elat-tgtlat)“2 + (elong - tgtlong)2)/pirad*60; 

range_error = [range_error,error]; 

end %%% for i= 1: 360 

figure(1) 

plot(flongout/pirad, flatout/pirad,’b*',elongmat,elatmat,'r+’,... 
tgtlongmat, flatout/pirad,'y-',index,tgtlatmat,'y-') 

xlabel(‘Longitude (deg) East’), ylabel('Latitude (deg)') 


axis([0,3,-1.5,1.5]) 
print -deps fig4ch 
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figure(2) 


plot(range_error) 
title(‘distance error (Nmi)’) 
print -deps bz3fig2 


figure(3) 

plot(elatmat,'+') 

title(‘Predicted Target Latitude (deg) Target Located at 0 deg') 
hold 

plot(tgtlatmat) 

axis({0,1,-5,5]) 


figure(4) 
plot(elongmat,'+') . 
title(Predicted Target Longitude (deg) Target Located at 1.5 deg’) 
hold 
plot(tgtlongmat) 
axis([0,1,-5,5]) 


figure(5) 

plot(itout) 

title(Iterations Required to Make Estimate’) 
axis([0,i1,0,max(itout)+2]) 
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APPENDIX B. COBURN KALMAN FILTER CODE 


This appendix is a representative case of Coburn’s FORTRAN code written in 


MATLAB. 

ToT To %0Q% CKFAVG.M ToT % To % 
To% 

%o 

% Mgj S.P. Jones, USMC 

% This function isa MATLAB adaptation of the FORTRAN 

% program developed by L. Laddie Coburn in June, 1972. 

%o flong is aircraft position longitude 

%o flat is aircraft position latitude 

Yo vel is the aircraft velocity 

%o azz is the raw doa azimuth bearing measurement 

% wt is the weighting factor for the measurement(not used) 
% 1 is the index denoting the number of the present measurement 
% ncuts is the total number of measurements 

% toa is time of arrival 

% elong is the estimated emitter longitude in radians 

% elat is the estimated emitter latitude in radians 

% acc is df accuracy 1 sigma usually 10 degrees 

clear all 


%%% Initiate the program 


randn(‘seed',2) | % set seed generator 

ncuts = 150; 

index = [1:ncuts-1]; 

degperrad = 180/pi; % conversion factor from deg-->rad 


elatmatc = zeros(1,ncuts-1); 
elongmatc = zeros(1 ,ncuts-1); 
rngerrsumc = zeros(1,ncuts- 1); 
total = 10; 


for loops = 1:total 


flatout = []; 
flongout = []; 
azzout = []; 


elatoutc = []; 
elongoutc = []; 
mgerroroutc = []; 
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flong! = 0/degperrad; % 0 degrees long (rad) 
flatl = -1.5/degperrad; % 1.5 deg S. (rad) 

flong = flong1; 

flat = flatl; 


tgtlong = 1.5/degperrad; % 1.5 deg E. (rad) 

tgtlat = 0/degperrad; % at the equator 
tgtlatout = tgtlat/degperrad*ones(1,ncuts-1); 
tgtlongout =1.5*ones(1,ncuts-1); . 


acc = 10/degperrad; % df accuracy is about +/- 10 degrees 

vel = 360; % AIC velocity is 360 knots 

rngest = 1.414*90; % initial range estimate always 150 
toa = 0; 


fori=I:ncuts %travel 180 NMi measuring every 12 sec 
dferror = acc*(randn); 
azz = (pi/2)-atan((tgtlat-flat)/(tgtlong-flong)) + dferror; 


azzout = [azzout,azz]; 
azzdeg = azz*degperrad; 


ifi== 1 


p11 = 10000; % P matrix 
pi2=0; 

p21 =0; 

p22 = 10000; 

ql1 =0; % Q matrix 

ql2 =0; 

q22 = 0; 

oldd11 = 1; % D 
oldd12 = 0; 

oldd21 = 0; 

oldd22 = 1; 

wli=1; 

wl2=0; 

w21 =0; 

w22=1; 

R = (acc)2; % (10 degree std dev)42 = var 


gl =p11/(p11+R); % G 
g2 = p12/(p11+R); 


thtd = azzdeg; % thtd will be the filtered df cut 

thtd1 = azzdeg; % thtd1 will the smoothed first cut 

thetadot = vel/rngest*sin(azz)*degperrad/3600; % theta dot deg/sec 
thetadot1 = thetadot; % smoothed first theta dot 

toa = 0; | 

timel = toa; % dummy variable for computing elpased time 
t=0; % elapsed time between cuts 
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else 


t = toa - time]; 


timel = toa; 
tt = t/1000; 


gill =tt; 
q12 = tt4(1.5); 


q22 = tt’2; 





Zoreset time! for next cut 


%o%% KALMAN RECURSION EQUATIONS %%%% 


newp11 = p11*(1-g1) + 2*p12*t - (pl2*g1 + p11*g2)*t + ... 


(p22 - p12*g2)*t*2 + q11; 
p 


newp12 = p12*(1-g1) + (p22 - p12*g2)*t + q12; 


newp22 = p22 - p12*g2 + q22; 


gl =newp11/(newp11 + R); 
g2 = newp12/(newp11 + R); 


tptd = thtd + thetadot*t; 


e = azzdeg - tptd; 


thtd = tptd + g1*e; 


thetadot = thetadot + g2*e; 


% filtered df cut 


% bearing rate 


Jo%oTo SMOOTHING EQUATIONS FOR FIRST CUT 2%%% 


%% Coburn relied on inverse 


%% matrix fumction that is not supported by CMS-2. 


P= [p11 
p21 


Pinv = inv(P); 


p12 
p22); 


pini1 = Pinv(1,1); 
pin12 = Pinv(1,2); 
pin22 = Pinv(2,2); 


gp11 = q11*pin11 + q12*pin12; 
qp12 = qi1*pin12 + q12*pin22; 
qp21 = q12*pin11 + q22*pin12; 
qp22 = q12*pin12 + q22*pin22; 


dil = oldd11*(1-qp11) - oldd12*qp21 + oldd11*qp21*t; 
d12 = oldd12*(1-qp22) - oldd11*qp12 + oldd11*(qp22-1)*t; 
d21 = oldd21*(1-qp11) - oldd22*qp21 + oldd21*qp21*t; 
d22 = oldd22*(1-qp22) - oldd21*qp12 + oldd21*(qp22-1)*t; 
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oldd11 =d11; 
oldd12 = d12; 
oldd21 = d21; 
oldd22 = d22; 


thtd1 = thtdl + (d11*g1 + d12*g2)*e; 
thetadot] = thetadot] + (d21*g1 + d22*g2)*e; 


pli =newp11; 
p12 = newp12; 
p22 = newp22; 


thti = thtd/degperrad; 
tht1i = thtd1/degperrad; 


tat = tan(tht11) - tan (thti); 

wav = (flatl + flat)/2; 

tila = ((flong-flong1)*cos(wav) + flat1*tan(tht1i) - flat*tan(thti))/tat; 
tla = ((flong-flong1)*cos(tila) + flat]*tan(tht1i) - flat*tan(thti))/tat; 


elat = tla*degperrad; 


tlo = flong + (tla - flat)*tan(thti)/cos(tila); 
elong = tlo*degperrad; 


elatoutc = [elatoutc,elat]; 
elongoutc = [elongoutc,elong]; 


mgerror = sqrt((tgtlat-tla)*2 + (tgtlong-tlo)“2)*degperrad*60; 
rngerroroutc = [rngerroroutc,rngerror]; 


end %(if iz=1 


oldlat = flat; 
oldlong = flong; 
toa = toa + 12; 
flat = flat + 12/60*6/60/degperrad; 
flatout = [flatout,flat]; 


end % fori= 1:91 

elatmatc = elatmatc + elatoutc; 
elongmatc = elongmatc + elongoutc; 
rngerrsumc = rngerrsumc + rngerroroutc; 
end % for loops = 1: 

elatavgc = elatmatc/total; 

elongavgc = elongmatc/total; 

rmgerravgc = rngerrsumc/total; 

figure(1) 


plot(index,elatavgc,'*b',index,tgtlatout,'r’) 
title(‘Estimated Latitude’) 
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axis([O,ncuts,-5,5]) 


figure(2) 
plot(index,elongavgce,'*b',index,tgtlongout,'r') 
title‘Estimated Longitude’) 
axis([O,ncuts,0,3]) 


figure(3) 
plot(mgerravgc) 
title(Error (Nmi)') 
axis([O,ncuts,0,100]) 
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APPENDIX C. MILLS’ KALMAN FILTER CODE 


This appendix is a representative case of Mills’ FORTRAN code written in 


MATLAB. 


To To Vo Vo No MKFAVG.M ToT %%Q% 

ToQTo TN 
% % 
% Maj S.P. Jones, USMC 

% This function is a MATLAB adaptation of the FORTRAN 

% program developed by E.H. Mills in March, 1973. 


%o flong is aircraft position longitude 

Jo flat is aircraft position latitude 

%o vel is the aircraft velocity 

To azz is the raw doa azimuth bearing measurement 
% i denotes the number of the present measurement 
Jo toa is time of arrival 

% elong is the estimated emitter longitude in radians 
Jo elat is the estimated emitter latitude in radians 

% acc is df accuracy 1 sigma usually 10 degrees 


%  thtis angle theta True Azimuth - True Heading I have 
Yo modeled the system around an aircraft heading 000 and 


%o eliminated this calculation 

% tht] is the smoothed first azimuth 

To thetadot is the bearing rate in rad/sec 

% thetadot] is the smoothed first bearing rate in rad/sec 
clear all 


%%% Initiate the program 


randn(‘seed',2) % set seed generator 

ncuts = 150; 

index = [1:ncuts-1]; 

degperrad = 180/pi; % conversion factor from deg-->rad 


elatmat = zeros(1,ncuts-1); 
elongmat = zeros(1,ncuts-1); 
rngerrsum = zeros(1,ncuts-1); 


total = 1; 


for loops = 1:total 


loops 

flatout = []; 
flongout = []; 
azzout = []; 
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elatout = [}; 
elongout = []; 
rmgerrorout = [{]; 


flong1 = 0/degperrad; % 0 degrees long (rad) 
flatl =-1.5/degperrad; % 1.5 deg S. (rad) 

flong = flong1; 

flat = flat1; 


tgtlong = 1.5/degperrad; % 1.5 deg E. (rad) 

tgtlat = O/degperrad; % at the equator 
tgtlatout = tgtlat/degperrad*ones(1,ncuts-1); 
tgtlongout =1.5*ones(1,ncuts-1); 


acc = 10/degperrad; % df accuracy is about +/- 10 degrees 
vel = 360; % AIC velocity is 360 knots 
rngest = 1.414*90; % initial range estimate always 150 


toa = 0; 

fori=Ii:ncuts  %travel 180 NMi measuring every 12 sec 
dferror = acc*(randn); 

azz = (pi/2)-atan((tgtlat-flat)/(tgtlong-flong)) + dferror; 


azzout = [azzout,azz]; 
azzdeg = azz*degperrad; 


ifi== 


pil = 10000; % P matrix 
p12 =0; 

p22 = 10000; 

ql1 =0; Yo Q matrix 

qg12 =0; 

q22 = 0; 

dil =1; % D 

wll=1; 

wl12=0; 

w21 = 0; 

w22 = 1; 

R = (acc)42; % (10 degree std dev)42 = var 


g1=pl11/(p11+R); % G 
g2 = p12/(p11+R); 


tht = azz; 
tht] = azz; 


thetadot = vel/rngest*sin(azz)/3600; ~ % theta dot deg/sec 
thetadot] = thetadot; % smoothed first theta dot 

toa = 0; 

time! = toa; % dummy variable for computing elpased time 
t=0; % elapsed time between cuts 
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else 
t= toa - timel; 
time! = toa; %reset time1 for next cut 
tt = t/1000; ) 
qil = tt; 
ql2 = tt4(1.5); 
q22 = tt2; 
%%% KALMAN RECURSION EQUATIONS %%%% 


newpl1 = p11*(1-g1) + 2*p12*t - (pl12*gi + pl1*g2)*t + ... 
(p22 - p12*g2)*t*2 + q11; 


newp12 = p12*(1-g1) + (p22 - p12*g2)*t + q12: 
newp22 = p22 - p12*g2 + q22; 


newgl = newp11/(newp11 + R); 
newg2 = newp12/(newp11 + R); 


tpt = tht + thetadot*t; 


e = azz - tpt; 
tht = tpt + g1*e; % filtered df cut 
thetadot = thetadot + g2*e; _ % bearing rate 


%%% SMOOTHING EQUATIONS FOR FIRST CUT 2%%% 
%% These equations are extracted from Edward H. Mills thesis 
%% that builds upon Coburn's work. Coburn relied on inverse 
%o% matrix fumction that is not supported by CMS-2, and so I 
%7% have used this group of equations as developed by Mills 
%% to better simulate operation onboard an actual aircraft. 


smooth! = 1 - p11*(1-g1)/R; 
smooth2 = smooth] *t; 
smooth3 = -p12 * (1-g1)/R; 
smooth4 = 1 + smooth3*t; 


neww11 = wll*smooth!l + w12*smooth?2; 
neww12 = w11*smooth3 + w12*smooth4; 
neww21 = w21*smooth1 + w22*smooth?; 
neww22 = w21*smooth3 + w22*smooth4; 


wil =newwll; 
wl12 = newwl12; 
w21 = neww21; 
w22 = neww22; 


d11 =d11 - (w1142)*(p11 + R)/(R412); 


thtl = tht] + (w11/R)*e; 
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thetadot] = thetadot1 + (w21/R)*e; 


pll =newp11; 
pi2 = newp12; 
p22 = newp22; 


gl =newsg!l; 
g2 = newg2; 


thti = tht; 
thtli = thtl; 


tat = tan(tht1i) - tan (thti); 

wav = (flat] + flat)/2; 

tila = ((flong-flong1)*cos(wav) + flat1*tan(tht1i) - flat*tan(thti))/tat; 
tla = ((flong-flong1)*cos(tila) + flat] *tan(tht1i) - i la 
tlo = flong + (tla - flat)*tan(thti)/cos(tila); 


elat = tla*degperrad; 
elong = tlo*degperrad; 


elatout = [elatout,elat]; 
elongout = [elongout,elong]; 


mgerror = (sqrt((tgtlat-tla)*2 + (tgtlong-tlo)*2))*degperrad*60; 
rngerrorout = [rngerrorout,mgerror]; 


end =. %(if i== 


oldlat = flat; 

oldlong = flong; 

toa = toa + 12; 

flat = flat + 12/60*6/60/degperrad; 
flatout = [flatout,flat]; 


end % fori= 1:ncuts 


elatmat = elatmat + elatout; 
elongmat = elongmat + elongout; 
rngerrsum = mmgerrsum + rmngerrorout; 


end % for loops = 1: 


elatavg = elatmat/total; 
elongavg = elongmat/total; 
mngerravg = mgerrsum/total; 


figure(1) 

plot(index,elatavg,'*b’,index,tgtlatout,'r') 

Zotitle(‘Estimated Latitude 150 Cuts 10+/-Degree Accuracy Mills Kalman Filter’) 
xlabel(‘Cuts'), ylabel(‘Latitude (Deg)') 

axis({0,ncuts,-5,5]) 

Yoprint -deps mkf12f1 
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figure(2) 

plot(index,elongavg,'*b',index,tgtlongout,'r') 

Zotitle(‘Estimated Longitude 150 Cuts 10+/-Degree Accuracy Mills Kalman Filter’) 
xlabel('Cuts'), ylabel(Longitude (Deg)') 

axis([O,ncuts,0,3]) 

Yoprint -deps mkf12f2 


figure(3) 

plot(index,rngerravg,'~',index tmgerrorout,'*') 

legend(‘1000X Simulated AVG','1X Simulation’) 

Ztitle(‘Error (NMi) 150 Cuts 10+/-Degree Accuracy Mills Kalman Filter’) 
xlabel(‘Cuts’), ylabel(‘Error (NMi)') | 

axis([{O,ncuts,0,200]) 

Yoprint -deps mkf12f3 


figure(4) 

plot(elongout,elatout,'*g') 

ylabel('Latitude (Deg)'), xlabel(‘Longitude (Deg)’') 

grid 

Zotitle(‘Mills Kalman Filter Predictions 150 Cuts +/- 10 Degree Accuracy’) 
axis((0,3,-.5,1.5]) 

Zoprint -deps mkf12f4 


figure(5) 

plot(index,rngerravg,'-',index,rngerrorout,'*') 

legend('1000X Simulated AVG','1X Simulation’) 

Zotitle(‘Error (NMi) 150 Cuts 10+/-Degree Accuracy Mills Kalman Filter’) 
axis({0O,ncuts,0,60]) 

xlabel(‘Cuts'),ylabel(Error (NMi)') 

Yoprint -deps mkf12f5 
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